﻿#include "gaulen.h"

// #include "p_complex.cpp"  // �����Ҫ�Ļ�����ȡ��ע��

#define EPS 3.0e-11

void gaulen(double x1, double x2, double* x, double* w, int n) {
    int m, j, i;
    double z1, z, xm, xl, pp, p3, p2, p1;

    m = (n + 1) / 2;
    xm = 0.5 * (x2 + x1);
    xl = 0.5 * (x2 - x1);

    for (i = 1; i <= m; i++) {
        z = cos(3.141592654 * (i - 0.25) / (n + 0.5));

        do {
            p1 = 1.0;
            p2 = 0.0;
            for (j = 1; j <= n; j++) {
                p3 = p2;
                p2 = p1;
                p1 = ((2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / static_cast<double>(j);
            }

            pp = n * (z * p1 - p2) / (z * z - 1.0);
            z1 = z;
            z = z1 - p1 / pp;
        } while (fabs(z - z1) > EPS);

        x[i] = xm - xl * z;
        x[n + 1 - i] = xm + xl * z;
        w[i] = 2.0 * xl / ((1.0 - z * z) * pp * pp);
        w[n + 1 - i] = w[i];
    }
}

/*complex ff(double x)
{
    return(exp(-x*x));
}

void main()
{
    double *x, *w;
    complex f = 0;
    int i, n = 1000, key;
    x = new double[n + 1];
    w = new double[n + 1];
    gaulen(0., 4.3, x, w, n);

    for (int i = 1; i <= n; i++) {
        std::cout << "\nx[" << i << "]=" << x[i] << "; w[" << i << "]=" << w[i];
    }

    std::cin >> key;

    for (i = 1; i <= n; i++) {
        f += w[i] * ff(x[i]);
    }

    std::cout << "f=" << f;
    std::cin >> key;

    delete[] x;
    delete[] w;
}*/
